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Introduction and Methodology 


Over 80 per cent of the silicon used by the electronics industry is 
grown from the melt by the Czochralski (CZ) crystal growth process, which 
is depicted schematically in Fig. 1 . A single crystal boule is slowly pulled 
from a pool of melt maintained by heating the outside of the crucible. 
Surface tension acts against the force of gravity to form a meniscus connect- 
ing the the growing crystal to the surface of the melt pool. The precise 
temperature gradient needed to sustain the solidification of nearly perfect 
crystalline solid is maintained by a cooler ambient above the crucible. 
The crystal and crucible are rotated to minimize thermal asymmetry in the 
system. In addition, as the melt level drops during growth, the crucible 
is raised by movement of the pedestal so that the solidification front remains 
in a specified region of the heater. 

Capillarity and heat transfer determine the shape and stability of crys- 
tals grown by the Czochralski method. We describe the solution of the first 
quasi-steady, thermal-capillary model of the CZ growth process which directly 
predicts the dependence of the crystal radius, melt-solid interface, and 
melt/gas meniscus on growth parameters and thermal boundary conditions in 
the furnace. The equations governing heat transfer are presented in Fig. 2., 
and the expressions for determining the radius and interface location are 
shown in Fig. 3. We model heat transfer in the system as conduction-dominated 
and include latent heat generation at the solidification interface. Radiation 
and convection accounts for heat exchange between the melt and crystal and 
the thermal surroundings. The thermal boundary conditions used in the calcu- 
lations presented here represent an idealized, but qualitatively realistic 
thermal environment for CZ crystal growth. The heating of the melt through 
the crucible is modeled by specifying the temperature of the melt adjacent 
to the crucible wall as a constant without coupling this value to the actual 
power input to the crucible from the surrounding heater. Heat flows out 
of the melt through the bottom of the crucible, according to a specified 
heat transfer coefficient set to model energy loss to the pedestal and the 
surroundings. By itself, this heat transfer system represents a slightly 
nonlinear conductionradiation problem and its solution would be straight- 
forward except that the shape of the melt and crystal are unknown and are 
aprt of the solution. 

The locations of the interfaces are coupled to the heat transfer in 
the system by setting the shapes of the regions for conduction and by determin- 
ing the surface areas for radiative and convective exchange with the ambient. 
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The melt/solid interface is determined by the location of the melting point 
isotherm. The shape of the melt meniscus is found from the Young-Laplace 
equation of capillary statics, the condition for three-phase contact at the 
melt/solid/gas tri-junction, the wetting condition at the crucible wall, 
and by the specification of the instantaneous volume of the melt. The radius 
of the crystal is determined so that the three-phase (melt/crystal/gas) equil- 
ibrium contact angle for growth is satisfied (Surek and Chalmers 1975). 

These relations form a complete set of coupled partial differential 
equations for this complex nonlinear free-boundary problem. Only Ettouney 
et al. (1983) have solved a thermal-capillary model with a similar degree 
of complexity. As detailed in Fig. 4., We solve the field and boundary equat- 
ions simultaneously by a finite-element/Newton scheme for the temperature 
field in melt and crystal, the shapes of the meniscus and melt-solid inter- 
face, and the radius of the crystal. Besides giving quadratic convergence 
in each of these unknowns, the Newton method is the basis for 
computer-implemented analyses of parametric sensitivity and stability 
(Yamaguchi et al. 1984). The details of the finite element algorithm are 
presented in a recent manuscript (Derby et al. 1985) for the modeling of 
liquid encapsulated Czochralski growth of gallium arsenide. 

Results 

Results for the thermophysical properties characteristic of CZ Silicon 
growth in a small-scale system are shown in Fig. 5. along with a specified 
ambient temperature distribution. On the left representative isotherms are 
shown for the melt and crystal. Note that heat flows into the melt from 
the crucible wall and exits the melt both upward through the surface of the 
melt and through the crystal and downward out of the bottom of the crucible. 
The shape of the meniscus is shown along with the melt/solid interface which 
is concave into the crystal. For this set of parameters, the resultant crys- 
tal radius was 3.8 cm. On the right, the ambient temperature distribution 
for this case is shown. The almost step decrease in the surrounding tempera- 
ture field that occurs above the top of the hot crucible was modeled by the 
ambient temeperature profile shown in Fig. 5. The hot crucible has an increas- 
ingly important effect on the crystal size as the melt level drops and the 
outer crystal surface directly views this distribution. Similar ambient 
temperature distributions were used for the decreasing melt volume calculat- 
ions of Figs. 10., 11., and 12. 

The response of the system to changes in the steady-state growth rate 
is detailed in Fig. 6. as the pull rate, represented by the dimensionless 
Peclet number, is varied from 0-20 cm/hr (0£Pe£0.4). All other parameters 
were held constant and the ambient distribution was taken to be uniform at 
1262K. The crystal radius decreased with increasing pull rate because of 
the increased latent heat release at the interface. This effect has been 
observed experimentally. The release of latent heat also was responsible 
for the increasing curvature of the solidification interface. At higher 
growth rates the interface became increasingly concave into the crystal. 

Some computational aspects of these calculations are presented in Figs. 
7 and 8. The converged finite-element meshes for each of the calculations 
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are shown in Fig. 8. Meshes with 76 element were used in all calculations 
presented in this report. The mesh follows the shape and position of the 
changing interfaces for each case. Since the mesh was based on the location 
of the free-boundaries, it converged as part of the entire solution with 
every iteration. The iteration histories for these runs is detailed in Fig. 
8. The convergence was quadratic, as expected for full Newton's method; 
see Ettouney and Brown (1983). 

Calculations are presented in Fig 9. for different equilibrium growth 
angles for the wetting of the crystal by the melt. The angle is varied from 
the reported 11 degree value for silicon to M5 and 90 degrees. There were 
two major effects: the crystal radius decreased with increasing angle, and 
the shape of the melt-solid interface changed from concave into the crystal 
to convex into the melt as the angle increased from 11 to 90 degrees. The 
last case for an equilibrium angle of 90 degrees corresponded to the assumpt- 
ion used by most previous investigators of a flat melt surface joining the 
crystal at right angles. This simplification obviously lead to erroneous 
results. Accurate modeling of heat flow in this region necessitates inclus- 
ion of the melt meniscus and the proper equilibrium angle. 

The next series of figures represents the interaction of the temperature 
field and interfaces with the melt volume. Several effects are present as 
the melt level drops. The amount of heat entering the melt from the crucible 
sides decreased because the surface area in contact with the wall decrea- 
sed. Also, the heat loss out of the bottom of the crucible was increasingly 
important in the shallower melt, and the hot crucible wall was exposed. 
The exposed wall allowed for radiative heat transfer from the hot crucible 
wall directly into the growing crystal. First, we examine only the effect 
of the dropping melt level without the complicating effects of radiative 
heat transfer. The case of decreasing melt volume with an axially uniform 
ambient temperature of 1262 K is shown in Fig. 10. Since the amount of heat 
entering the melt decreased as the melt level drops, the radius of the crystal 
increased. The effect of the heat loss through the bottom of the crucible 
was manifested in a change of the interface shape from concave to sigmoidal. 
This transition is known as interface flipping and commonly has been attri- 
buted to transitions in the natural convection within the melt. Since our 
calculations are for a conduction-dominated system, interface flipping was 
caused solely by conductive heat losses through the melt and the changing 
shape of the melt pool. Another related effect of the crucible heat loss 
was the formation of a section of undercooled melt at the bottom of the cru- 
cible for the shallow melt depth (V m »0.5ir). In practice, this region would 
be a small cap of solid silicon at the bottom of the crucible which can grow 
large enough with decreasing melt depth to touch the upper melt/solid inter- 
face. This is occasionally observed in practice, and is known as "crystal 
bumping". 

When we add radiation from the exposed wall in the form of the specified 
axial temperature distribution discussed previously, we get the situations 
depicted in Figs. 11. and 12. Even though the heat flow into the melt 
decreased with the dropping melt, the radiative transfer of heat from the 
exposed wall into the crystal can prevent the growth in the radius of the 
crystal as the melt level dropped. This effect is seen for the case with 
melt volume V m «o.5ir for the shallow crucible in Fig. 11. and for both cases 
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for the deep crucible in Fig. 12. Radiation from the crucible wall to the 
crystal is so overpowering for the deeper crucible case of Fig. 12 that a 
non-zero crystal radius did not exist for the V m =o.5ir case. Radiative inter- 
change between wall and crystal also magnified the interface flipping effect. 
With heat entering the side of the crystal, the melt/solid interface was 
pushed even further into the melt. The solid cap at the bottom of the cru- 
cible was present for the melt volume V m =0.5it shown in Fig. 11. 

The changes in radius for decreasing melt volume are shown together 
in Fig. 13. For the uniform ambient, the crystal radius increased steadily 
with decreasing melt volume. When an ambient temperature distribution was 
introduced to include the crucible wall the increasing effect of radiation 
from the exposed wall the radius show either a maximum with decreasing volume 
or decreased monotonically, depending on the depth of the crucible. For 
the shallow crucible, the crystal radius initially increased as the melt 
dropped, but decreased as more of the crucible wall was exposed to the view 
of the crystal. This effect is even more dramatic for the deep crucible: the 
radius decreased from the start and rapidly melted away so that no steady- 
state growth was possible after some minimum volume was reached. This last 
effect is often seen in practice where the crystal becomes very difficult 
to grow as the melt level drops below some critical value. 

Summary 

We have demonstrated the success of efficiently calculating the tempera- 
ture field, crystal radius, melt mensicus, and melt/solid interface in the 
Czochralski crystal growth system by full finite-element solution of the 
governing thrmal-capillary model. The model predicts realistic response 
to changes in pull rate, melt volume, and the thermal field. The experi- 
mentally observed phenomena of interface flipping, bumping, and the difficulty 
maintaining steady-state growth as the melt depth decreases are explained 
by model results. These calculations will form the basis for the first quant- 
itative pitcure of CZ crystal gowth. The accurate depiction of the melt 
meniscus is important in calculating the crystal radius and solidification 
interface. The sensitivity of our results to the equilibrium growth angle 
place doubt on less sophisticated attempts to model the process without inclus- 
ion of a the meniscus. 

Quantitative comparison with experiments should be possible once more 
representation of the radiation and view factors in the thermal system and 
the crucible are included. Extensions of our model in these directions are 
underway. 
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Noraencl ature 


Bij, Biot number - ^Rq/^s for surface j; 

Bo, Bond number - gR c 2 p/o; 

Cp, heat capacity [J/kgK]; 

£2, unit axial vector; 

g, gravitational constant [m/s 2 ]; 

hj, heat transfer coefficient for surface j [W/m 2 K]; 

height of interface i from bottom of crucible [m] ; 
hj, dimensionless interface heights « hj/R c ; 
k, thermal conductivity [W/mK]; 

Ki, temperature- dependent thermal conductivity ratio - ki/k s for region i; 
ni, normal vector to surface i; 

Pe, Peclet number * VpR 0 p s Cp s /k s ; 

r, radial coordinate measued from center of crucible [m]; 
r, dimensionless radial coordinate « r/R c ; 

R, crystal radius [m]; 

R, dimensionless crystal radius = R/R c ; 

R c , crucible radius [m]; 

Rai, Radiation number = 0* iRc^f^^^s ^ or surface i; 

S, Stefan number = Hf/Cp s Tf; 

T, temperature [K]; 

Tf, solidification temperature [K]; 

V m , volume of melt [m 3 ]; 

V m , dimensionless volume of melt « V m /R c 3 ; 

Vp, steady-state crystal pull rate [m/s]; 

z, axial coordinate measured from bottom of crucible [m]; 

z, dimensionless axial coordinate « z/R c ; 
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Greek symbols 

Hf, heat of fusion [J/kg-mole]; 
p, density difference across melt meniscus; 
j, emissivity of surface j; 

li, dimensionless reference pressure » Pi.daturr/ Pi8 R c5 
4>, equilibrium wetting angle; 

р, density [kg/m 2 ]; 

o, surface tension [ J/m] ; 

o*, Stef an-Boltzmann constant [W/m 2 K^]; 

0 , dimensionless temperature * T/Tf; 

Subscripts 
a, ambient; 

с, crucible; 

i s j, numerical indeces denoting surface, domain, or equation; 
ra, melt; 

s, solid. 
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Czochralski Crystal Growth 
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Steady-State Thermal-Capillary Model 
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CAPILLARITY AND INTERFACES 
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FINITE ELEMENT/NEWTON METHOD 


■ Piacewiee-continuous polynomials uead to 
approx imato temperature field (2d) and 
interface shapes (Id) 


s Adaptive mash defined by free boundaries 


s Galerkin method and elemental isoparametric 
mopping used to generate eat of nonlinear 
algebraic equations 


■ Newton's method used to solve aquation set. 
giving quadratic convergence to solution 


ISOPARAMETRIC MAPPING 



-a 
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Parameters— Silicon 




Ambient Temperature 


206 



Finite-Element Meshes for Different Pull Rates 
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Decreasing Melt Volume — Shallow Crucible 
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Decreasing Melt Volume — Deep Crucible 
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Radius as a Function of Melt Volume 
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DISCUSSION 


ELWELL: Is it too complicated to plug convection in the Czochralski system? 

BROWN: It depends on your goal. We are working on convection, and we are 

working on convection in other systems that are quite a bit simpler 
than this. To come up with a good algorithm for doing convection in 
a large Czochralski system means you have to handle convection that 
is in the never-never land between laminar and turbulent flow, which 
is essentially chaotic. It is too complicated at present to plug 
that kind of convection model into a calculation like this and hope 
to get back out some kind of a control strategy. Those are 
supercomputer-type calculations, unless you are willing to put in a 
turbulent model, which we are not willing to do right now. 

DUDUKOVIC: You use a static value of the contact angle in all of your 

calculations, which is of course fine because presumably the pulling 
rates are so low. What do you do about the value of that angle in 
presence of the Marangoni effect? 

BROWN: If you really believe that is a thermal physical property of what is 

right down at the intersection between melt/solid and gas then you 
would believe that the bulk fluid mechanics should not have a lot of 
effect on them. That concept has been shown to have been true not 
in solidification systems but in normal wetting conditions as long 
as the contact line is not moving at high rates of speed. 

DUDUKOVIC: At the moment that’s the best assumption one can make, but isn't 

there still some doubt about it under dynamic conditions? 

BROWN: If it's a microscopic property of the physics of the tri-junction I 

don’t think there is any scaling you can do to show that a bulk flow 
can affect it. For example, gravity cannot affect it. The same 
kind of argument you would use to say that the surface tension is 
independent of gravity is used and, in fact it’s one order of 
magnitude worse. 

KIM: In your aluminum segregation study, what was the effective segregation 

coefficient that you have measured? 

BROWN: One. 

KIM: What was the solid boundary layer thickness? 

BROWN: There is no boundary layer thickness because in the EFG system that is 
a big problem. The flow field in the die is well developed and 
there is no effective segregation in the system. Any aluminum that 
comes up through that die has to go into the crystal. It has no 
place else to go. 

KALEJS: It is quite true for K e ff equal to unity that there are a lot of 

impurities that ended up very close to the meniscus, but this is 
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true for a lot of the shaped crystal growth systems like web. We 
have looked at this and found that the change in surface tension 
with impurity concentration goes opposite to the temperature 
difference, so you may expect some kind of cancellation due to 
impurity effects, particularly in K e ff equal to unity where a lot 
of impurities do end up in the meniscus. 

SUREK: Regarding Professor Dudukovic's question: the contact angle was 

determined over many orders of magnitude of growth rate, essentially 
from zero growth rate to very fast growth rate without much 
variation, so it seems to be a property of the silicon at a triple 
junction. My question is if you showed points or regions where 
there is no existence of a steady-state solution, would that be the 
limit of stability, in your opinion? 

BROWN: The EFG calculations give you the growth-rate limit for steady-state 

solutions for growing a ribbon with that die at that K e ff. We had 
not found similar limits in the Czochralski system, but I think all 
one needs to do is play around with the thermal ambient to actually 
find growth-rate limits in that system also. I don't think people 
worry about it as much because they concentrate more on high-quality 
material and less on growth speed. 

SUREK: In what system was it that you found no solution? 

BROWN: In the Czochralski system. Essentially that’s associated with just 

excessive heat transfer from the ambient directly into the crystal, 
which essentially just blocks up the heat so that you can’t get it 
out of the melt. We are not sure whether or not that is essentially 
just the crystal growing to zero diameter or whether or not there 
really is a limit point there. We haven't refined that enough to 
find out. 

WITT: How do you handle the complication that normally is encountered when you 

grow dislocation-free silicon? In conventional growth, you end up 
with the peripheral interface facet, which places the critical 
portion of the crystal into a supercooled state, and we know that 
you lose, for one reason or another, the dislocation-free growth 
configuration. The instantaneous consequence is a significant 
change in crystal diameter, whibh means that effectively the growth 
does not occur in critical portions, particularly at the periphery, 
which is a critical point way away from the melting point. My 
second question: a very important element is the rate of crucible 
consumption or the melt/crucible interaction or the mass transfer 
related to oxygen and impurities. Are you looking at these factors? 

BROWN: We have not been. Obviously, if you take the same segregation model we 

can dissolve the crucible with it if you want. We have not been 
doing it. In answer to your first question, interface faceting is a 
fantastic piece of physics that could be put in here if I knew the 
orientation and undercooling law for doing it. There's nothing that 
keeps us from putting in a faceted growth mode at the interface. It 
would be very interesting to do, and there are some mathematical 
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reasons why it's very nice to do in the finite-element formulation. 
You could very easily come up with an interface that has facets 
along part of it and equilibrium growth on other parts. If someone 
would stick their neck out and give me the undercooling law I'd love 
to do it, because I think it would be very interesting. I think 
; that introduces an additional non-linearity, which I think can give 
you catastrophic behavior, which is what you referred to with the 
loss of the facet in the growth of the crystal. Do you know of 
someone that has written down this nucleation law? 

WITT: It' is still, very controversial in details. The theory is worked out. 
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